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SUMMARY 

Clipping of narrow extrema and distortion of smooth profiles is a well-known 
problem associated with so-called "high-resolution” nonoscillatory convection 
schemes. In this report, a strategy is presented for accurately simulating highly 
convective flows containing discontinuities such as density fronts or shock waves, 
without distorting smooth profiles or clipping narrow local extrema. The convection 
algorithm is based on non-artificially-diffusive third-order upwinding in smooth 
regions, with automatic adaptive stencil expansion to (in principle, arbitrarily) 
higher order upwinding locally, in regions of rapidly changing gradients. This is 
highly cost-effective because the wider stencil is used only where needed - in isolated 
narrow regions. A recently developed universal limiter assures sharp monotonic 
resolution of discontinuities without introducing artificial diffusion or numerical 
compression. An adaptive discriminator is constructed to distinguish between 
spurious overshoots and physical peaks; this automatically relaxes the limiter near 
local turning points, thereby avoiding loss of resolution in narrow extrema. 
Examples are given for one-dimensional pure convection of scalar profiles at 
constant velocity. 
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The University of Akron, Akron, Ohio 44325. 

♦♦Doctoral student, College of Engineering, The University of Akron, Akron, Ohio 44325. 


INTRODUCTION 


Accurate simulation of highly convective flows continues to be one of the most 
challenging problems in computational mechanics. The inability to adequately 
simulate simple scalar profiles even in the superficially straight-forward case of one- 
dimensional pure convection at constant velocity has been referred to as the 
"ultimate embarrassment” for computational mechanics [1], Classical (second-order 
central) schemes are able to handle very smooth profiles fairly satisfactorily; i.e., 
when the convected variable has small curvature (spatial second derivatives) - or, 
more specifically, when the shortest Fourier components have wavelengths of many 
mesh-widths. But a sudden change in gradient (involving locally high curvature or 
short wavelengths) excites spurious unphysical oscillations due to inherent spatial 
(third-derivative) dispersion terms in the truncation error. By contrast, first-order 
upwinding [2] can resolve a step profile monotonically, giving a spreading error- 
function in the idealized case; but the inherent spatial second-derivative (artificial 
diffusion) terms in the truncation error so totally overwhelm any physical diffusion 
that even smoothly varying profiles are unrecognizable after a short time. 

Second-order upwinding [3] represents a significant improvement. In this case, 
leading order truncation error is again a dispersive third derivative but, compared 
with second-order central, fourth- (and higher even-) derivative dissipation terms 
may be stronger under certain conditions, thus tending to somewhat dampen disper- 
sive oscillations. By averaging second-order upwinding with the canonical explicit 
second-order central scheme [4], Fromm [3] proposed a convection algorithm in 
which the leading phase error of the upwind scheme approximately cancels the 
lagging phase error of the central scheme. However, using the same computational 
stencil, it is a simple matter to eliminate the third-derivative dispersion term 
entirely. 

The QUICKEST scheme [5] is the canonical explicit third-order upwind con- 
vection algorithm of this type; physical diffusion terms are modelled to a consistent 
order. In this case, the leading truncation-error term is a dissipative (but not diffu- 
sive) spatial fourth derivative. This is an important attribute. ^Numerical experi- 
ence over the past decade or more has led to the following rule of thumb: leading 
truncation error in a convection algorithm should be dissipative (an even spatial 
derivative) rather than dispersive (odd derivative) but should be of higher order than 
other modelled physical terms such as diffusion (especially if these terms are sup- 
posed to be small). Thus, even-order schemes fail the first criterion (central schemes 
are worse in this respect because dissipative terms are small or non-existent). First- 
order upwinding fails the second criterion. Third-order upwinding is therefore the 
lowest order method satisfying both criteria for the convection-diffusion equation. 
This is especially important for convection-dominated flows, as will be seen in the 
model test problems shown later. 

Unfortunately, third-order upwinding does not seem to have been widely 
adopted by the computational-fluid-dynamics (CFD) community. Apparently for no 
other reasons than tradition and personal inertia, many CFD researchers regard 
second-order central schemes as the norm. This is appropriate for physical problems 
dominated by even spatial derivatives, such as diffusion, wave-motion, and elas- 
ticity, but is singularly inappropriate for the /irsf-derivative convection term in fluid 
mechanics [6]. A similar situation occurs with higher order odd derivatives such as 
the spatial third derivative in the Korteweg-de Vries equation [7]. Another large 
segment of the CFD community (especially in heat-transfer applications) has 
adopted first-order upwinding as the "safest” form of convection modelling. This 
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apparently stems from the influence of the Imperial College school through descen- 
dants of the TEACH code [8] using the "Hybrid” convection scheme [9] or variants of 
the so-called "exponential difference scheme” [10], including the "power-law’ 
approximation described in a well-known textbook [11]. The difficulty with such an 
approach is very simple; using these methods, one is not simulating the correct 
physical problem! A very-high-Reynolds-number (convection-dominated) problem is 
arbitrarily replaced by an unphysical problem in which the effective component grid 
Reynolds number can never be greater than 2, whereas the physical grid Reynolds 
number might be hundreds or thousands (or effectively infinite). 

There is also a serious logical flaw in such methods: sophisticated (and expen- 
sive) multi-equation turbulence models are usually used at each time or iteration 
step simply as a diagnostic to switch off their own effects in the turbulent-transport 
terms of the governing equations, replacing these terms with artificial viscosity or 
diffusivity [12]. Specifically, the turbulence model is used to calculate the (turbu- 
lent) component grid Reynolds (or Peclet) number at each control-volume face; when 
this quantity exceeds 2 (in the case of Hybrid) or about 6 (for the exponential or 
power-law schemes), physical turbulent transport is replaced by artificial numerical 
viscosity (or diffusivity). The short-comings of first-order methods are well known, 
having been shown theoretically [13-15] and through comparative benchmark 
problems [16-17]. Their continued wide-spread popularity is perhaps another 
manifestation of the ultimate embarrassment. 

Although third-order upwinding is a more rational basis for CFD, certain 
fundamental problems remain. Perhaps the most obvious is the fact that unphysical 
overshoots and undershoots are excited each side of a step discontinuity in purely 
convective flow [5]. And, of course, short-wavelength resolution is limited by the 
third-order accuracy. Switching to higher order (upwind) schemes can increase 
accuracy of resolution but cannot reduce the overshoot problem (in fact, it gets 
slightly worse). Monotonic high resolution schemes can be constructed for simu- 
lating step discontinuities. This is usually achieved by using a subtle nonlinear 
blending of a (traditional) second-order central base scheme with first-order 
upwinding (introducing enough positive artificial diffusion locally to maintain 
monotonicity) and first-order downwinding (with inherent negative artificial diffu- 
sion to locally enhance artificial "compression”). This strategy is the basis of so- 
called "shock-capturing” and "total-vanation-diminishing (TVD)” schemes [18], and 
can be related to earlier multi-step flux-correction first/second-order schemes [19] 
and single-step nonlinear flux-limiting schemes such as the second-order MUSCL 
scheme [20] - a monotonized version of Fromm’s method - and the third-order 
EULER scheme based on exponential upwinding [21,22]. Unfortunately, the best 
step-resolution (supercompressive) schemes such as Superbee [23], Super-C, or 
Hyper-C [24], tend to convert smooth profiles into a series of steps and plateaus. In 
other words, the negative diffusivity responsible for artificial compression of discon- 
tinuities tends to concentrate moderate-curvature regions into localized unphysical 
sharp changes in gradient. All such first/second-order schemes strongly clip 
narrow extrema. 

Recently, a universal limiter (UL) has been designed which can be applied to 
arbitrarily high order transient interpolation modelling (TIM) of the advective 
transport equations (ATE). This "ULTIMATE” CFD scheme [24] has several attrac- 
tive properties. Step resolution is monotonic and can be made competitive with the 
best supercompressive shock-capturing schemes without introducing artificial dif- 
fusivity or viscosity, simply by locally increasing the order of accuracy, using adap- 
tive stencil expansion (so that the location of the wider stencil automatically moves 
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along with high-curvature regions of the profile). Since numerical compression is 
avoided, smooth profiles are not corrupted. However, one annoying problem 
remains: very narrow extrema are slightly clipped relative to what can be achieved 
with an unlimited higher order (upwind) scheme. Near sudden changes in gradient, 
the limiter needs to be "on”, but near local physical extrema (but not spurious over- 
shoots), it needs to be switched "off” automatically. The problem is one of pattern 
recognition. This paper describes a simple adaptive discriminator which can identify 
well-defined local narrow physical extrema and automatically switch off the 
universal limiter in such regions, but keep it activated in regions where unphysical 
overshoots or very short- wavelength spurious oscillations would otherwise occur. 

The overall strategy is as follows: 

(i) In "smooth” regions (identified by small values of local absolute curvature of 
the converted variable), the unlimited third-order upwind QUICKEST scheme 
is used; this accounts for the overwhelming bulk of the flow domain - especially 
in multidimensional flows. 

(ii) In relatively large-gradient or strong-curvature regions, automatic adaptive 
stencil expansion occurs locally, using (in principle, arbitrarily) higher order 
upwinding with the universal limiter activated. 

(iii) Near well-defined local extrema, identified by the automatic discriminator, the 
limiter is switched off, thus allowing an appropriate degree of resolution 
depending on the narrowness of the extremum. 

The next section summarizes the basic ideas underlying the universal limiter. 
Then four challenging scalar test profiles are considered under conditions of pure 
one-dimensional convection at constant velocity. Results for a number of solution 
methods are shown in order to demonstrate some of the difficulties mentioned above 
with respect to various manifestations of the ultimate embarrassment. Construction 
of the automatic discriminator is then briefly described. Finally, results are shown 
for a cost-effective variable-order scheme consisting of QUICKEST in smooth 
regions, ULTIMATE seventh-order upwinding in moderate-curvature regions and 
ULTIMATE ninth-order upwinding in high-curvature or high-gradient regions. 


UNIVERSAL LIMITER 

Details of the ULTIMATE scheme can be found elsewhere [24]; however, a brief 
description of the universal limiter is given here for convenience. Figure 1 shows a 
one-dimensional control volume, suggesting the local behaviour of the convected 
variable <|> in terms of locally normalized variables [22] 

4 > ( 1 ) 

where is the (unnormalized) upstream node value and <J> W the downstream 
value. Note from Equation (1) and the figure that, in terms of normalized variables, 

$»t, = 0 and $ u = 1 (2) 

Let <Jy be^ the normalized face value on the downstream control- volume face; 
similarly, <j> u is the normalized upstream face value. Figure 1 shows interpolative 
constraints on $y and §i u ;i.e., 
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(3) 


$c n = 1 


and 


0 ^ $ u ^ $ c " (4) 

where <|> c n is the normalized central node value. The inequalities given by (3) are 
two necessary conditions on <$y prescribed by the universal limiter when local 
behaviour is monotonic across tne three nodes shown in Figure 1, i.e., the two nodes 
straddling the face in question and the adjacent upstream node. However, this is not 
enough to guarantee computational monotonicity. For this, consider the explicit 
update step for §» c under conditions of constant velocity (positive to the right in 
Figure 1). In terms of normalized variables, 

$ c n+1 =$c n (5) 

where c = uAf/ Ax is the Courant number. 


In order to maintain monotonicity (locally), <J> C must satisfy 

^ $ c " +1 ^ ( 6 ) 

Taking conservative worst-case conditions implies 

0 = V +1 ^ 1 ' : ’(7) 

The right-hand inequality is assured by the interpolative constraints on and <J> U ; 
the left-hand inequality implies, using Equation (5), 

$ f ^ + $ c "/c (8) 

Once again, a^worst-case estimate for ( = 0) results in an additional (time-step) 
constraint on 

$ f %$ c "/c (9) 

Inequalities (3) and (9) constitute the universal limiter constraints on with 
respect to when 4> c n is within the monotonic range: 

0 S 1 (10) 

Outside of this range (i.e., <|> c n < 0 or v> 1), various strategies can be devised. 
Figure 2 shows one simple possibility: 

$ = $ c n (nonmonotonic range) (11) 

together with the monotonicity constraints. This is a graphical portrayal of the 
universal limiter. Using this limiter, the procedure is as follows: 

(i) Construct an explicit estimate for the unnormalized control-volume face value, 

by any (in principle, arbitrarily high order) method. 

(ii) Compute the corresponding normalized value, <$y, according to Equation (1). 
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(iii) If the point (4> c n , $/) liesjvithin the shaded region of Figure 2, proceed with the 
unadjusted face value, 

(iv) If the point lies outside this region, replace §y with the nearest allowable §y 
value at the same 

(v) Reconstruct $ f = + <$y (<J> 0 - 

(vi) Repeat for each control-volume face. 

Note that this strategy is significantly different from that used in currently 
popular shock-capturing and TVD schemes [18]. In such schemes, using the present 
normalized- variable terminology, 

$ f = c$ c * + (l-c)§y" (12) 

where is a single-valued monotonic function of <$L n passing through the points 
(0,0), (0.5, 0.75), and (1,1). This restricts such methods to second (or, at best) third 
order accuracy in space, and second order in time. Figure 3 shows the functional 
form of §y n (<l>£ n ) for three well-known methods: Roe’s Minmod and Superbee [23], 
and van Leers MUSCL [20], As commonly used in shock-capturing schemes, the 
TVD limiter [25] can be described by 

^ min(2$' c n , 1) (13) 

in the monotonic range, with $y n = elsewhere. These constraints, together with 
Equation (12), represent much more restrictive conditions than those of the uni- 
versal limiter [24]. An important point to stress is that the universal limiter can be 
used with arbitrarily high order (in both space and time) explicit methods whereas 
the usual TVD strategy is restricted to essentially second order. 


CONVECTION OK SCALAR PROFILES 

Consider one-dimensional constant-velocity pure convection (i.e., no diffusion) 
of a scalar profile initially given by 

<J>(x,0) = <J> o (x) (14) 

After time t , the exact solution is the initial profile translated by a distance ut, 
where u is the convecting velocity; i.e., 

4>(x,f) = $ o (x-ut) (15) 

Figure 4 shows the exact solution of four test profiles under uniform convection to 
the right. Discontinuity resolution is tested by a simple step rather than the square 
wave used by some authors [25] (at best, a square wave simply gives twice as much 
information as a step; for the more dispersive schemes, oscillations produced by the 
square-wave’s step-up can interfere with those produced by the step-down, thus 
producing a confusing interaction pattern). Smooth behaviour is represented by an 
isolated sine-squared wave 20Ax wide; there is no discontinuity in value or 
gradient, but there is a discontinuity in curvature at each side of the profile. A semi- 
ellipse 20Ax wide is very challenging because of a combination of sudden and 
gradual changes in gradient. A narrow Gaussian with o= 1.94A* has been chosen 
to test narrow peak resolution. 
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In the following simulations, the Courant number for all cases is 


(16) 


c = 0.45 

Calculations are run for 100 time-steps, so that exact solutions should have shifted 
45 mesh-widths to the right. The exact solutions are shown for reference in each 
case. All of the schemes have been programmed in conservative control-volume 
form, so that "mass” is conserved (to machine accuracy) in all cases. Upstream 
(<|> = l) and downstream (<J> = 0) boundaries are sufficiently far away (beyond the 
region shown in the figures) that they do not interfere with the solutions. To get 
some idea of how far the profiles have moved, note that the three maxima are each 
30Ax apart. Thus, for example, the initial position of the Gaussian peak was half- 
way between the final positions of the sine-squared and semi-ellipse maxima. 

Figure 5 shows the results of using first-order upwinding. Because of strong 
artificial diffusion (spatial second derivative) in the truncation error, the computed 
profiles have strongly interacted. Individually, the step simulation would look like a 
spreading error-function, whereas the other profiles soon degenerate into spreading 
(and amplitude-decaying) Gaussians, since all short-wavelength components in the 
original profiles are quicxly damped out. 

If one defines the grid Peclet number as 

P A = uAx/D (17) 

where D is the physical diffusion coefficient, it is clear that, the physical problem 
under consideration corresponds to infinite P A . As is well known [26], first-order 
upwinding for infinite P A is indistinguishable from any higher order method with 
an effective grid Peclet number given by 

P A * = 2/(1 -c) (18) 

In other words, for transient problems, first-order upwinding introduces artificial 
numerical diffusion of the form 

D nwn = uAx(l-c)/2 (19) 

Such artificially diffusive results should be recognized as being totally unacceptable. 
It is indeed rather surprising that such inaccurate methods are still the industry 
standard in many branches of computational engineering, especially numerical 
convective heat transfer [27,28]. 

Figure 6 shows results for the explicit second-order central Lax-Wendroff [4] 
scheme. In this case, there is a large phase lag due to the spatial third-derivative 
term in the truncation error. Only the smooth (sine-squared) profile bears any 
resemblance to the corresponding exact solution. Clearly, for highly convective 
flows, the explicit Lax-Wendroff method is very disappointing for all but the 
smoothest of profiles. The situation is not improved by using explicit second-order 
upwinding, shown in Figure 7. In this case, the third-derivative truncation-error 
terms are of opposite sign to those of the Lax-Wendroff method, resulting in phase- 
lead dispersion. For Courant numbers near 0.5, the individual second-order upwind 
results are almost a reflection of the Lax-Wendroff profiles. 
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Implicit second-order convection schemes are not necessarily better than 
explicit methods. The so-called "Crank-Nicolson-type” convection scheme [26] - 
denned by analogy with the well-known method for tne diffusion equation [29] - is 
actually worse than explicit methods in the purely convective case, because there are 
strong phase-lag dispersion terms and no (even-derivative) dissipation terms in the 
truncation error. The highly oscillatory nature of this method is seen in Figure 8. 
The second-order implicit linear -finite-element method [26], shown in Figure 9, 
again has no dissipation but gives somewhat better results than other second-order 
methods because the leading third-derivative dispersion term is smaller. 

Clearly, it would seem somewhat more felicitous to design a high-convection 
scheme in which the oscillatory third-derivative dispersion term in the truncation 
error were entirely eliminated. This was achieved over ten years ago in a simple 
explicit formulation known as the QUICKEST scheme [5]. In this case, leading trun- 
cation error is a small (fourth-derivative) dissipation term which strongly damps any 
dispersive tendencies of the non-zero fifth-derivative term without corrupting the 
modelled physics, viz.: the absence of physical (second-derivative) diffusion terms. 
Figure 10 shows a dramatic improvement over the dispersive second-order methods. 
Note particularly the good phase behaviour (which is relatively insensitive to 
Courant number). For reference, Figure 1 1 shows results for an explicit fourth-order 
central method [24]. Phase-lag dispersion in this case stems from the leading fifth- 
derivative truncation error (which is only lightly damped by higher order dissipa- 
tion). Results are only slightly better than those of the best second-order method. 

At the risk of belabouring the point, the above results seem to indicate fairly 
clearly that third-order upwinding represents a much more rational basis for CFD 
[6] than first- or second-order methods or higher order central schemes. As 
mentioned previously, there are still two serious short-comings, both evident in 
Figure 10: (i) unphysical undershoots or overshoots near regions of rapid change in 
gradient (high curvature); (ii) a limit to short- wavelength resolution. The first of 
these deficiencies can be overcome by applying the universal limiter; Figure 12 
shows the resulting ULTIMATE QUICKEST computation. Step resolution is 
monotonic and comparable to that of MUSCL, Figure 13, a limited form of Fromm’s 
method. (One should note that at c = 0.5, Fromm’s method is equivalent to 
QUICKEST [24] ). But it is clear that MUSCL tends to clip and flatten local extrema 
more strongly. This is directly linked to the more restrictive conditions imposed by 
the commonly used "TVD” limiter, given by (13). 

The second deficiency of third-order upwinding can, of course, be corrected by 
using higher order (upwind) methods. The strategy proposed here is to use an 
appropriate higher order method locally, with third-order upwinding in relatively 
"smooth” (low-curvature) regions. For reference, Figure 14 shows results for 
(unlimited) ninth-order upwinding [24] used globally; ninth-order appears to be 
necessary to fully capture the peak of the narrow Gaussian. But, of course, this is an 
unlimited scheme, and undershoots and overshoots are excited near discontinuities. 
ULTIMATE ninth-order, Figure 15, eliminates the oscillations and gives very tight 
step-resolution, but again introduces slight clipping at extrema. If one could 
somehow arrange to automatically switch off the limiter near local extrema (while 
keeping it active near discontinuities), a very desirable convection scheme would 
result. This is the job of an adaptive discriminator which can distinguish between 
well-defined local physical extrema and short-wavelength spurious oscillations. 
Note that long- wavelength numerical oscillations, such as those occurring in second- 
order simulations - which might otherwise "confuse” the discriminator - present no 
problem (because they are excluded from the basic scheme). 
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CONSTRUCTION OF THE DISCRIMINATOR 


The following algorithm attempts to distinguish between artificial numerical 
peaks, such as those potentially occurring near the step or the base of the semi- 
ellipse, and true physical extrema, such as those of the exact profiles in the test cases 
studied. Artificial peaks would be associated with short-wavelength numerical 
oscillations, with rapidly changing value, gradient, and curvature; as mentioned 
above, long-wavelength numerical oscillations are excluded from the third/higher 
order algorithm. If a local extremum is associated with short-wavelength oscilla- 
tions, the discriminator chooses the limited algorithm; however, if the curvature is of 
the same sign in adjacent regions, the discriminator relaxes the limiter constraints 
at the extremum and the immediate upstream and downstream adjacent nodes. 

In setting up the convective fluxes (for the right-face of each CV cell), consider 
an increasing-i "DO-loop” sweep. For the discriminator, choose a stencil of seven 
points: d>.,, <b „, <J>, , <F, 4> 1+9 , Then compute the differences between 

each pair of consecutive points: 


D1 = (<!>< . 2 , - . 3 ), D2 = ($. , - 4> j 2 ) D6 = ($ I+3 - $; + 2 ) 


( 20 ) 


For convenience, assume there is a local maximum; a minimum requires reversal of 
some of the subsequent inequalities. Limiter constraints are active unless otherwise 
stated. The algorithm proceeds as follows:. 

(i) Check if D1 and D2 are both positive and D3 and D4 both negative; if not, 
go to step (iii); if they are, then: 

(ii) Check if |D2| < |D1| and |D3| < |D4| ; if true, switch to an unlimited version at 
the current i-value and skip the remaining steps; otherwise: 

(iii) Check if D2 and D3 are both positive and D4 and D5 both negative; if not, 
go to step (v); if they are, then: 

(iv) Check if |D3| < |D2| and |D4| < |D5| ; if true, switch to an unlimited version at 
the current i-vaiue and skip the remaining steps, otherwise: 

(v) Check if D3 and D4 are both positive and D5 and D6 both negative; if not, 
keep limiter constraints active; if they are, then: 

(vi) Check if |D4| < ID3| and ]D5| < |D6| ; if true, switch to an unlimited version at 
the current i-vaiue; if not, proceed with the limiter constraints active. 

Figure 16(a) shows a sketch of a case in which the discriminator would keep the 
universal limiter activated; in Figure 16(b), the limiter constraints would be 
removed at each of the points shown by a hollow circle. The discriminator routine is 
by-passed unless one of $ + { , or 4>j . is a local extremum. As noted above, limiter 
relaxation occurs (if at all) in groups of three points, with the discrete extremum in 
the middle. 


ADAPTIVE STENCIL EXPANSION 


The ultimate convection scheme proposed here uses the unlimited QUICKEST 
(third-order upwind) algorithm in smooth non-steep regions; i.e., wherever the 
average (right-face) absolute "curvature” 

CURVAV = 0.5 I*,,, - 4>, +1 • + $,.,1 (2D 
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and "gradient” 

GRAD = |4> 1+I - 4>J (22) 

are both less than pre-assigned thresholds, and provided [24] 

0.2 ^ I'c ^ 0.8 (23) 

If CURVAV exceeds THC1 or the larger THC2 , the algorithm automatically 
switches locally (at the CV face in question) to ULTIMATE seventh- or ninth-order 
upwinding, respectively. If GRAD exceeds THG it locally switches directly to 
ULTIMATE ninth-order upwinding. This order-switching strategy is sketched in 
Figure 17. Also, if the limits expressed in (23) are locally exceeded, the limiter 
constraints are activated for the basic third-order scheme. Near local extrema, the 
discriminator automatically relaxes the limiter constraints, as described in the 
previous section. Clearly, the (dimensional) threshold constants need to be chosen 
carefully so as to capture the desired degree of accuracy without "overkill” in terms 
of cost; this requires some experimentation for optimization; a more universal 
(nondimensional) procedure is currently being explored. 

Adaptive stencil expansion is a very cost-effective strategy because the wider- 
stencil computations are performed only where needed - at a relatively small 
number of grid points in narrow regions, by definition. This is even more effective in 
two and three dimensions. Figure 18(a) shows results for the complete algorithm; 
the small arrows show where the automatic discriminator has relaxed the limiter 
constraints (at this particular time-step). Figure 18(b) shows the corresponding 
distribution of GRAD , and 18(c) that of CURVAV. Finally, Figure 18(d) shows the 
local order of the algorithm to be used (at the next time-step) for the threshold 
constants chosen. It should be clear that the locations of the limiter-relaxation 
points and the order-switches move along with the profiles automatically, as time 
progresses. 


CONCLUSION 

A cost-effective strategy for high-convection modelling has been introduced. 
The algorithm is based on the (third-order upwind) QUICKEST scheme in smooth 
regions, as this is the lowest order method in which the leading truncation error is 
dissipative but not diffusive. QUICKEST’s phase error is relatively low, but 
overshoots or undershoots can be excited by the unlimited scheme near discon- 
tinuities, and short- wavelength resolution would, of course, be limited to third order. 
A more sophisticated scheme is therefore devised in which monotonicity is 
guaranteed by a universal limiter which can be applied to any order of accuracy. 
Appropriate degrees of higher order resolution are automatically introduced locally, 
using adaptive stencil expansion controlled by monitoring local (absolute) gradient 
and curvature of the convected variable. Potential clipping of narrow extrema is 
avoided by automatic relaxation of the limiter constraints in local regions based on 
pattern-recognition decisions of an adaptive discriminator. By judicious choice of 
threshold constants, the overall method can produce extremely accurate results on a 
coarse mesh at relatively little additional cost above that of the base third-order 
scheme. The same principles can be applied to multidimensional flows and non- 
linear systems. Work is proceeding on these applications. 
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FIGURE 1 


FIGURE 2 



Locally monotonic behaviour across a control-volume cell in terms 
of normalized variables. 



Universal limiter constraints in the normalized-variable diagram. 
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FIGURE 3 Normalized variable diagrams for Superbee (heavy lines), Minmod 

(dashed), and MUSCL. 
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FIGURE 6 


Lax-Wendroff results compared with exact solution. 
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FIGURE 7 Second-order upwind results compared with exact solution. 
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FIGURE 9 Linear-finite-element results compared with exact solution. 



FIGURE 10 Third-order upwind (QUICKEST) results compared with exact 

solution. 
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FIGURE 1 1 Fourth-order central results compared with exact solution. 



FIGURE 12 ULTIMATE QUICKEST results compared with exact solution. 
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FIGURE 13 Results using MUSCL compared with exact solution. 



FIGURE 14 Unlimited ninth-order upwinding. 
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FIGURE 15 


ULTIMATE ninth-order upwinding. 



(a) (b) 

FIGURE 16 Action of the discriminator, (a) Limiter switched on. 

(b) Limiter switched off at three central points. 
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FIGURE 18 


Third/higher-order adaptive scheme with discriminator. 

(a) Computed results, (b) Distribution of GRAD relative to THG. 

(c) Distribution of CURVAV relative to THC1 and THC2. 

(d) Order of the algorithm to be used locally in next time-step. 
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